Geometry operations

Note

This is an R Notebook based on chapter 5 of the book Geocomputation with R written by Robin Lovelace, Jakub Nowosad, and Jannes Muenchow. Minor changes have been done in order to produce outputs from each chunk. Example datasets have also been changed.

Prerequisites

This session uses the same packages as the previous session on spatial-operations:

library(sf)
library(raster)
library(tidyverse)
library(spData)
library(spDataLarge)

Introduction

The previous G4D sessions have demonstrated how geographic datasets are structured in R and how to manipulate them based on their non-geographic attributes and spatial properties. This session extends these skills.

After reading it — and attempting the exercises at the end — you should understand and have control over the geometry column in sf objects and the geographic location of pixels represented in rasters.

Geometric operations on vector data

This section is about operations that in some way change the geometry of vector (sf) objects. It is more advanced than the spatial data operations presented in the previous session because here we drill down into the geometry: the functions discussed in this section work on objects of class sfc in addition to objects of class sf.

Simplification

Simplification is a process for generalization of vector objects (lines and polygons) usually for use in smaller scale maps. Another reason for simplifying objects is to reduce the amount of memory, disk space and network bandwidth they consume: it may be wise to simplify complex geometries before publishing them as interactive maps. The sf package provides st_simplify(), which uses the GEOS implementation of the Douglas-Peucker algorithm to reduce the vertex count. st_simplify() uses the dTolerance to control the level of generalization in map units [see @douglas_algorithms_1973 for details].

Simplification is also applicable for polygons. This is illustrated using countries, representing borders of contiguous countries in the world. As we show in a previous section, GEOS assumes that the data is in a projected CRS and this could lead to unexpected results when using a geographic CRS. Therefore, the first step is to project the data into an apropriate projected CRS:

library(rnaturalearth)
sf_countries <- st_as_sf(ne_countries())
#
sf_sam <- sf_countries %>% filter(subregion == "South America")
tsf_sam = st_transform(sf_sam, 31991)
plot(st_geometry(tsf_sam), col = sf.colors(12, categorical = TRUE), border = 'black', 
     axes = TRUE)

st_simplify() works equally well with projected polygons:

sam_s1 = st_simplify(tsf_sam, dTolerance = 100000)  # 100 km
plot(st_geometry(sam_s1), col = sf.colors(12, categorical = TRUE), border = 'black', 
     axes = TRUE)

A limitation with st_simplify() is that it simplifies objects on a per-geometry basis. This means the ‘topology’ is lost, resulting in overlapping and ‘holy’ areal units illustrated. ms_simplify() from rmapshaper provides an alternative that overcomes this issue.

# proportion of points to retain (0-1; default 0.05)
library(rmapshaper)
sam_s2 = rmapshaper::ms_simplify(tsf_sam, keep = 0.10,
                                          keep_shapes = TRUE)

Finally, the visual comparison of the original dataset and the two simplified versions shows differences between the Douglas-Peucker (st_simplify) and Visvalingam (ms_simplify) algorithm outputs:

Polygon simplification in action, comparing the original geometry of South America countries with simplified versions, generated with functions from sf (center) and rmapshaper (right) packages.

Polygon simplification in action, comparing the original geometry of South America countries with simplified versions, generated with functions from sf (center) and rmapshaper (right) packages.

Centroids

Centroid operations identify the center of geographic objects. Like statistical measures of central tendency (including mean and median definitions of ‘average’), there are many ways to define the geographic center of an object. All of them create single point representations of more complex vector objects.

The most commonly used centroid operation is the geographic centroid. This type of centroid operation (often referred to as ‘the centroid’) represents the center of mass in a spatial object. Geographic centroids have many uses, for example to create a simple point representation of complex geometries, or to estimate distances between polygons. They can be calculated with the sf function st_centroid() as demonstrated in the code below, which generates the geographic centroids of south american countries to the River Amazonas.

# rivers
rivers <- ne_download(scale = 110, type = 'rivers_lake_centerlines', category = 'physical')
## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/l_/1nv5l0dx78g_d3k3_t85lrjr0000gn/T/RtmpRSmoG2", layer: "ne_110m_rivers_lake_centerlines"
## with 13 features
## It has 31 fields
## Integer64 fields read as strings:  scalerank ne_id
sf_rivers <- st_as_sf(rivers)
amazonas <- sf_rivers %>% filter(name == "Amazonas")
tamazonas = st_transform(amazonas, 31991)
plot(st_geometry(tsf_sam), col = 'gray', border = 'black', 
     axes = TRUE)
plot(st_geometry(tamazonas), col = "blue", add = TRUE)

sam_cent = st_centroid(tsf_sam)
unique(sam_cent$name)
##  [1] "Argentina"    "Bolivia"      "Brazil"       "Chile"       
##  [5] "Colombia"     "Ecuador"      "Falkland Is." "Guyana"      
##  [9] "Peru"         "Paraguay"     "Suriname"     "Uruguay"     
## [13] "Venezuela"
plot(st_geometry(tsf_sam), col = 'gray', border = 'black', 
     axes = TRUE)
plot(st_geometry(tamazonas), col = "blue", add = TRUE)
plot(st_geometry(sam_cent), col = "red", cex=0.5, pch=21, add = TRUE)

dist2amz <- st_distance(sam_cent, tamazonas)
dist2amz
## Units: [m]
##              [,1]
##  [1,] 2355334.652
##  [2,]  789974.313
##  [3,]  907679.758
##  [4,] 2546746.120
##  [5,]  885047.217
##  [6,]  698225.417
##  [7,] 4298189.709
##  [8,]  829048.577
##  [9,]    4666.618
## [10,] 1688508.946
## [11,]  674563.487
## [12,] 2561146.947
## [13,] 1110104.805
peru_cent <- sam_cent %>% filter(name == "Peru")
peru2amz <- st_distance(peru_cent, tamazonas)
peru2amz
## Units: [m]
##          [,1]
## [1,] 4666.618

Other types of centroid exist, including the Chebyshev center and the visual center. We will not explore these here but it is possible to calculate them using R.

Buffers

Buffers are polygons representing the area within a given distance of a geometric feature: regardless of whether the input is a point, line or polygon, the output is polygon. Unlike simplification (which is often used for visualization and reducing file size) buffering tends to be used for geographic data analysis. How many points are within a given distance of this line? Which demographic groups are within travel distance of this new shop? These kinds of questions can be answered and visualized by creating buffers around the geographic entities of interest.

Before creating buffers, let’s create an sf spatial objects with countries’ capital cities in South America. First, let’s create an spatial object representing South America (SA) using the sf spatial union operator:

SA <- st_union(tsf_sam)
SA
## Geometry set for 1 feature 
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -3041628 ymin: -6302509 xmax: 2319657 ymax: 1467729
## epsg (SRID):    31991
## proj4string:    +proj=utm +zone=22 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## MULTIPOLYGON (((-688115.7 -5979162, -688115.7 -...

Now, let’s retrieve a rnaturalearth dataset with country capitals. Then, create the cities object using the spatial subsetting operation learned in the last week:

cities <- ne_download(scale = 110, type = 'populated_places', category = 'cultural')
## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/l_/1nv5l0dx78g_d3k3_t85lrjr0000gn/T/RtmpRSmoG2", layer: "ne_110m_populated_places"
## with 243 features
## It has 119 fields
## Integer64 fields read as strings:  wof_id ne_id
sf_cities <- st_as_sf(cities)
tsf_cities = st_transform(sf_cities, 31991)
sam_cities <- tsf_cities[SA, op = st_within]

Now, we are ready to create buffers of different sizes (50 and 500 km) surrounding the river Amazonas (and to find which cities are within each of them).

The command st_buffer() requires at least two arguments: an input geometry and a distance, provided in the units of the CRS (in this case meters):

amz_buff_50km = st_buffer(tamazonas, dist = 50000)
amz_buff_500km = st_buffer(tamazonas, dist = 500000)
Buffers around the River Amazon of 50km (left) and 500km (right). Note the colors, which reflect the fact that one buffer is created per geometry feature.

Buffers around the River Amazon of 50km (left) and 500km (right). Note the colors, which reflect the fact that one buffer is created per geometry feature.

Now, let’s find which cities are within a 500 km distance from the River Amazonas:

amz_cities <- sam_cities[amz_buff_500km, op = st_within]
amz_cities$NAME
## [1] "La Paz" "Lima"

Affine transformations

Affine transformation is any transformation that preserves lines and parallelism. However, angles or length are not necessarily preserved. Affine transformations include, among others, shifting (translation), scaling and rotation. Additionally, it is possible to use any combination of these. Affine transformations are an essential part of geocomputation.

For example, shifting is needed for labels placement, scaling is used in non-contiguous area cartograms, and many affine transformations are applied when reprojecting or improving the geometry that was created based on a distorted or wrongly projected map.

The sf package implements affine transformation for objects of classes sfg and sfc.

sam_sfc = st_geometry(tsf_sam)

Shifting moves every point by the same distance in map units. It could be done by adding a numerical vector to a vector object. For example, the code below shifts all y-coordinates by 10 kilometers to the north but leaves the x-coordinates untouched (left panel on the Fig. @ref(fig:affine-trans)).

sam_shift = sam_sfc + c(0, 1000000)

Scaling enlarges or shrinks objects by a factor. It can be applied either globally or locally. Global scaling increases or decreases all coordinates values in relation to the origin coordinates, while keeping all geometries topological relations intact. It can by done by subtraction or multiplication of asfg or sfc object.

sam_scale0 = sam_sfc * 0.5
plot(st_geometry(tsf_sam), col = 'gray', 
     axes = TRUE)
plot(st_geometry(sam_shift), col = "azure", add = TRUE)
plot(st_geometry(sam_scale0), col = "beige", add = TRUE)

Local scaling treats geometries independently and requires points around which geometries are going to be scaled, e.g. centroids. In the example below, each geometry is shrunk by a factor of two around the centroids (central panel on the Fig. @ref(fig:affine-trans)). To achieve that, each object is firstly shifted in a way that its center has coordinates of 0, 0 ((nz_sfc - nz_centroid_sfc)). Next, the sizes of the geometries are reduced by half (* 0.5). Finally, each object’s centroid is moved back to the input data coordinates (+ nz_centroid_sfc).

sam_centroid_sfc = st_centroid(sam_sfc)
sam_scale = (sam_sfc - sam_centroid_sfc) * 0.5 + sam_centroid_sfc
plot(st_geometry(tsf_sam), col = 'gray', 
     axes = TRUE)
plot(st_geometry(sam_scale), col = "beige", add = TRUE)
plot(st_geometry(sam_centroid_sfc), col = "azure", cex=0.5, pch=21, add = TRUE)

Rotation of two-dimensional coordinates requires a rotation matrix:

\[ R = \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \\ \end{bmatrix} \]

It rotates points in a counterclockwise direction. The rotation matrix can be implemented in R as: 

rotation = function(a){
  r = a * pi / 180 #degrees to radians
  matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
}

The rotation function accepts one argument a - a rotation angle in degrees. Rotation could be done around selected points, such as centroids (right panel on the Fig. @ref(fig:affine-trans)). See vignette("sf3") for more examples.

sam_rotate = (sam_sfc - sam_centroid_sfc) * rotation(30) + sam_centroid_sfc
Illustrations of affine transformations: shift, scale and rotate.

Illustrations of affine transformations: shift, scale and rotate. 

sam_scale_rotate = (sam_sfc - sam_centroid_sfc) * 0.25 * rotation(90) + sam_centroid_sfc

Finally, the newly created geometries can replace the old ones with the st_set_geometry() function:

sam_scale_sf = st_set_geometry(tsf_sam, sam_scale)
plot(st_geometry(sam_scale_sf), col = 'gray', border = 'black', 
     axes = TRUE)
plot(st_geometry(tamazonas), col = "blue", add = TRUE)

Clipping

Spatial clipping is a form of spatial subsetting that involves changes to the geometry columns of at least some of the affected features.

Clipping can only apply to features more complex than points: lines, polygons and their ‘multi’ equivalents. To illustrate the concept we will start with a simple example. Imagine you want to select not a 1000 km River Amazonas buffer nor a 500 km SA capitals buffer, but the space covered by both x and y. This can be done using the function st_intersection(), illustrated using objects named x and y which represent the two buffers:

x= st_buffer(tamazonas, dist = 1000000)
y = st_buffer(sam_cities, dist = 500000)

x_and_y = st_intersection(st_geometry(x), st_geometry(y))

plot(st_geometry(x), col = 'beige', border = 'black', 
     axes = TRUE)
plot(st_geometry(y), col = "azure", add = TRUE)
plot(x_and_y, col = "lightgrey", add=TRUE) # color intersecting area
plot(st_geometry(tamazonas), col = "blue", lwd=4, add=TRUE) # add river
plot(st_geometry(sam_cities), col = "red", cex=3, pch=20, add=TRUE) # add cities
Overlapping circles with a gray color indicating intersection between them.

Overlapping circles with a gray color indicating intersection between them.

Several times, It is relevant to obtain the difference between two geometries:

y_nor_x = st_difference(st_geometry(y), st_geometry(x))

plot(st_geometry(x), col = 'beige', border = 'black', 
     axes = TRUE)
plot(st_geometry(y), col = "azure", add = TRUE)
plot(y_nor_x, col = "salmon", add=TRUE) # color non intersecting area
plot(st_geometry(tamazonas), col = "blue", lwd=4, add=TRUE) # add river
plot(st_geometry(sam_cities), col = "red", cex=3, pch=20, add=TRUE) # add cities

Geometry unions

An spatial aggregation can silently dissolve the geometries of touching polygons in the same group. It is obtained using the st_union operator:

ins_cities <- sam_cities[y, op=st_within]
x_nor_y = st_difference(st_geometry(x), st_geometry(y))
cap_buf = st_buffer(ins_cities, dist = 500000)
un <- st_union(cap_buf, sam_s2,by_feature=TRUE)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
plot(st_geometry(un), col = "azure") # add union

It is also possible to get the __convex hull_ of a given geometry:

ins_cities <- sam_cities[y, op=st_within]
cvx <- st_convex_hull(st_union(ins_cities))
plot(st_geometry(sam_cities), col = "green", cex=2, pch=20) # add cities
plot(st_geometry(cvx), col = "grey", add=TRUE) # color intersecting area
plot(st_geometry(tamazonas), col = "blue", lwd=3, add=TRUE) # add river
plot(st_geometry(ins_cities), col = "red", cex=2, pch=20, add=TRUE) # add cities

What is going on in terms of the geometries? Behind the scenes, both aggregate() and summarize() combine the geometries and dissolve the boundaries between them using st_union().

The function can take two geometries and union them, as demonstrated in the code chunk below which creates a united western block incorporating Texas (challenge: reproduce and plot the result):

Geometric operations on raster data

In order to illustrate this topic, we start downloading a digital elevation model (DEM) directly from R. For such a purpose, we will use functionalities from the elevatr package. In order to download a DEM, we need to have an sp spatial object defining the area of interest.

The world dataset provided by rnaturalearth with the function ne_countries() is indeed an sp object:

sp_countries <- ne_countries()
sp_countries
## class       : SpatialPolygonsDataFrame 
## features    : 177 
## extent      : -180, 180, -90, 83.64513  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 63
## names       : scalerank,      featurecla, labelrank,  sovereignt, sov_a3, adm0_dif, level,              type,       admin, adm0_a3, geou_dif,     geounit, gu_a3, su_dif,     subunit, ... 
## min values  :         1, Admin-0 country,         2, Afghanistan,    AFG,        0,     2,           Country, Afghanistan,     AFG,        0, Afghanistan,   AFG,      0, Afghanistan, ... 
## max values  :         3, Admin-0 country,         7,    Zimbabwe,    ZWE,        1,     2, Sovereign country,    Zimbabwe,     ZWE,        0,    Zimbabwe,   ZWE,      1,    Zimbabwe, ...

Let’s select a region:

sp_deu_aus <- sp_countries[sp_countries$name == "Germany" | 
                         sp_countries$name == "Austria",]
sp_deu_aus
## class       : SpatialPolygonsDataFrame 
## features    : 2 
## extent      : 5.988658, 16.97967, 46.43182, 54.9831  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 63
## names       : scalerank,      featurecla, labelrank, sovereignt, sov_a3, adm0_dif, level,              type,   admin, adm0_a3, geou_dif, geounit, gu_a3, su_dif, subunit, ... 
## min values  :         1, Admin-0 country,         2,    Austria,    AUT,        0,     2, Sovereign country, Austria,     AUT,        0, Austria,   AUT,      0, Austria, ... 
## max values  :         1, Admin-0 country,         4,    Germany,    DEU,        0,     2, Sovereign country, Germany,     DEU,        0, Germany,   DEU,      0, Germany, ...

We will use the sp object to get our DEM using the get_elev_raster function. Note that we use zoom (z) = 2, that is, a very coarse resolution:

library(elevatr)
dem2 <- get_elev_raster(sp_deu_aus, z = 2, src = "aws")
dem2
## class       : RasterLayer 
## dimensions  : 514, 521, 267794  (nrow, ncol, ncell)
## resolution  : 0.176, 0.133  (x, y)
## extent      : -0.88, 90.816, -0.9177396, 67.44426  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : file5e5265dfb20 
## values      : -5365.321, 6044.875  (min, max)

Note that obtained DEM is only 634 rows by 516 columns. Note also that the cells are not square as they are usually.

raster::plot(dem2)

Geometric intersection

We can use the crop() function to “extract” a DEM raster object containing the cells whose midpoints overlap with the extent of a given sp spatial object

library(raster)
sp_deu <- sp_countries[sp_countries$name == "Germany",]
deu_dem <- crop(dem2, sp_deu)
deu_dem
## class       : RasterLayer 
## dimensions  : 57, 51, 2907  (nrow, ncol, ncell)
## resolution  : 0.176, 0.133  (x, y)
## extent      : 5.984, 14.96, 47.36126, 54.94226  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : file5e5265dfb20 
## values      : -63.142, 1811.796  (min, max)
plot(deu_dem)
plot(sp_deu, add = TRUE)

Let’s say that we need bioclimatic data for our country:

First, we need the centroid of our country:

library(sf)
sf_countries <- st_as_sf(sp_countries)
library(dplyr)
deu <- sf_countries %>% filter(name == "Germany")
deu_cent <- st_centroid(deu)
## Warning in st_centroid.sf(deu): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data
print(deu_cent)
## Simple feature collection with 1 feature and 63 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 10.28849 ymin: 51.13372 xmax: 10.28849 ymax: 51.13372
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   scalerank      featurecla labelrank sovereignt sov_a3 adm0_dif level
## 1         1 Admin-0 country         2    Germany    DEU        0     2
##                type   admin adm0_a3 geou_dif geounit gu_a3 su_dif subunit
## 1 Sovereign country Germany     DEU        0 Germany   DEU      0 Germany
##   su_a3 brk_diff    name name_long brk_a3 brk_name brk_group abbrev postal
## 1   DEU        0 Germany   Germany    DEU  Germany      <NA>   Ger.      D
##                     formal_en formal_fr note_adm0 note_brk name_sort
## 1 Federal Republic of Germany      <NA>      <NA>     <NA>   Germany
##   name_alt mapcolor7 mapcolor8 mapcolor9 mapcolor13  pop_est gdp_md_est
## 1     <NA>         2         5         5          1 82329758    2918000
##   pop_year lastcensus gdp_year                 economy
## 1       NA       2011       NA 1. Developed region: G7
##             income_grp wikipedia fips_10 iso_a2 iso_a3 iso_n3 un_a3 wb_a2
## 1 1. High income: OECD        NA    <NA>     DE    DEU    276   276    DE
##   wb_a3 woe_id adm0_a3_is adm0_a3_us adm0_a3_un adm0_a3_wb continent
## 1   DEU     NA        DEU        DEU         NA         NA    Europe
##   region_un      subregion             region_wb name_len long_len
## 1    Europe Western Europe Europe & Central Asia        7        7
##   abbrev_len tiny homepart                  geometry
## 1          4   NA        1 POINT (10.28849 51.13372)
tprec <- getData('worldclim', var='prec', res=0.5, lon=10.29, lat=51.13)
tprec
## class       : RasterStack 
## dimensions  : 3600, 3600, 12960000, 12  (nrow, ncol, ncell, nlayers)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 0, 30, 30, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names       : prec1_16, prec2_16, prec3_16, prec4_16, prec5_16, prec6_16, prec7_16, prec8_16, prec9_16, prec10_16, prec11_16, prec12_16 
## min values  :        2,        1,        1,        0,        0,        0,        0,        0,        0,         0,         1,         2 
## max values  :      291,      245,      271,      237,      244,      241,      229,      234,      273,       294,       266,       307

Note that we get a stack of 12 raster layers, each one representing monthly precipitation values (that is, an average value from a 30-year time series).

Let’s plot precipitation in November:

plot(tprec$prec11_16)
plot(sp_deu, add = TRUE)

Now, let’s crop the stack in order to get only the cells for our country:

deu_prec <- crop(tprec, sp_deu)
deu_prec
## class       : RasterBrick 
## dimensions  : 922, 1083, 998526, 12  (nrow, ncol, ncell, nlayers)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 5.991667, 15.01667, 47.3, 54.98333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : prec1_16, prec2_16, prec3_16, prec4_16, prec5_16, prec6_16, prec7_16, prec8_16, prec9_16, prec10_16, prec11_16, prec12_16 
## min values  :       21,       19,       26,       31,       42,       47,       52,       53,       37,        30,        29,        23 
## max values  :      169,      133,      148,      152,      151,      209,      191,      188,      143,       127,       162,       181

Note that, in the code chunk above, there is a vector-raster integration: we used the sp spatial object (vector data) as a “mask” to crop the precipitation raster layer.

plot(deu_prec$prec11_16)
plot(sp_deu, add = TRUE)

We have now two rasters:

par(mfrow=c(1,2))
plot(deu_prec$prec11_16, main="Precipitation - Nov -")
plot(sp_deu, add = TRUE)
plot(deu_dem, main = "DEM")
plot(sp_deu, add = TRUE)

We can notice that the DEM looks blocky because of its coarse spatial resolution.

Let’s review the spatial characteristics (i.e. dimension, resolution, extent,…) of the two raster we have got until now, the DEM and the precipitation:

deu_dem
## class       : RasterLayer 
## dimensions  : 57, 51, 2907  (nrow, ncol, ncell)
## resolution  : 0.176, 0.133  (x, y)
## extent      : 5.984, 14.96, 47.36126, 54.94226  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : file5e5265dfb20 
## values      : -63.142, 1811.796  (min, max)
deu_prec
## class       : RasterBrick 
## dimensions  : 922, 1083, 998526, 12  (nrow, ncol, ncell, nlayers)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 5.991667, 15.01667, 47.3, 54.98333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : prec1_16, prec2_16, prec3_16, prec4_16, prec5_16, prec6_16, prec7_16, prec8_16, prec9_16, prec10_16, prec11_16, prec12_16 
## min values  :       21,       19,       26,       31,       42,       47,       52,       53,       37,        30,        29,        23 
## max values  :      169,      133,      148,      152,      151,      209,      191,      188,      143,       127,       162,       181

We note that these raster datasets differ with regard to dimension, extent and spatial resolution.

Extent

Let’s start to solve this spatial mismatch using the extend function. It returns an Raster* object with a larger spatial extent.

ndeu_prec <- extend(deu_prec, deu_dem)
extent(ndeu_prec)
## class       : Extent 
## xmin        : 5.983333 
## xmax        : 15.01667 
## ymin        : 47.3 
## ymax        : 54.98333

Aggregation and disaggregation

In order to match spatial resolution we can either decrease (aggregate()) or increase (disaggregate) the spatial resolution of a raster. The raster package allows to do it on any of the two ways.

The aggregate function needs a fact (factor) parameter (i.e. 1, 2, 3, 4, 5) as well as a function parameter (i.e median, sum, max, min, …).

The dissagregate function also needs the factor parameter. In addition, we have to specify a method how to fill the new cells. There are two methods: nearest neighbor and bilinear. The first method gives all output cells the value of the nearest input cell, and hence duplicates values which leads to a pixelated output raster. The second method, in turn, is an interpolation technique that uses the four nearest cell centers of the input raster to compute an average weighted by distance.

While these functions are quite useful, our problem is more serious: the DEM cells are not square cells. After a search in the raster package documentation, I found the resample function which seems very useful here. Resample transfers values between non matching Raster* objects (in terms of origin and resolution).

ndeu_dem <- resample(deu_dem, ndeu_prec, method="bilinear")
ndeu_dem
## class       : RasterLayer 
## dimensions  : 922, 1084, 999448  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 5.983333, 15.01667, 47.3, 54.98333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : file5e5265dfb20 
## values      : -68.37913, 2018.066  (min, max)

Let’s review how the new DEM raster looks:

par(mfrow=c(1,2))
plot(ndeu_prec$prec11_16, main="Precipitation - Nov -")
plot(sp_deu, add = TRUE)
plot(ndeu_dem, main = "DEM")
plot(sp_deu, add = TRUE)

Raster extraction

Let’s imagine that we need to extract values of elevation at specific locations (i.e points):

deu_pts <- spsample(sp_deu, 20, type="regular")
deu_pts
## class       : SpatialPoints 
## features    : 21 
## extent      : 6.417043, 14.01875, 47.59302, 53.67439  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
plot(sp_deu)
plot(deu_pts, add = TRUE)

We can use the extract function to extract elevation values from ndeu_dem and assign the resulting values to a new column (elevation) in the deu_pts spatial object:

deu_pts$elevation <- raster::extract(ndeu_dem, deu_pts)
summary(deu_pts$elevation)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -4.40   57.59  243.40  287.58  435.86 1212.18

Rasterization

Sometimes, it is useful to convert a given spatial object from vector to raster (i.e. for spatial analysis of land suitability). The raster package provides such functionality. Here, we will explore just one case (of several possible cases as mentioned by Lovelace, Nowosad & Muenchow (2018)), the presence/absence raster.

raster_pts <- rasterize(deu_pts, ndeu_dem)
raster_pts
## class       : RasterBrick 
## dimensions  : 922, 1084, 999448, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 5.983333, 15.01667, 47.3, 54.98333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       :        ID, elevation 
## min values  :  1.000000, -4.400123 
## max values  :    21.000,  1212.181

Let’s plot the result:

par(mfrow=c(1,2))
plot(ndeu_dem, main = "DEM")
plot(sp_deu, add = TRUE)
plot(raster_pts$elevation, main="Points")
plot(sp_deu, add = TRUE)

We don’t notice, but the raster points are there.

nraster_pts <- aggregate(raster_pts, fact=50, fun= mean) 
nraster_pts
## class       : RasterBrick 
## dimensions  : 19, 22, 418, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 0.4166667, 0.4166667  (x, y)
## extent      : 5.983333, 15.15, 47.06667, 54.98333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       :        ID, elevation 
## min values  :  1.000000, -4.400123 
## max values  :    21.000,  1212.181

Let’s plot the aggregated raster points:

par(mfrow=c(1,2))
plot(ndeu_dem, main = "DEM")
plot(sp_deu, add = TRUE)
plot(nraster_pts$elevation, main="Points")
plot(sp_deu, add = TRUE)

Well, this is all for now.

Bis bald!!!